/*-------实验1--------*/
# include <stdio.h>
# include <math.h>
# define PI 3.14

void main()
{
	double Pan_d(double b, double a);//函数声明 

	int t; //循环变量位置
	double x_k[101], r_1[101], r_2[101], r_3[101], r_4[101], v_1[101], v_2[101], v_3[101], v_4[101];
	double b=200.0, l=400.0, x_o=1000.0, z_o=1000.0, a=90.0, i=90.0, M=2000.0, Q=2.67, z_k=0, san_g[101], san_x[101], san_z[101],G=6.67e-5, z_2, z_1; 
	
	FILE *fp; //文件指针
	fp = fopen("shujushuju.txt", "w+");
	/*----弧度转换----*/                                                                                                                                                                                                                                                                    
	a = (PI/180)*a;
	i = (PI/180)*i;
	/*------z1和z2的计算---------*/
	z_2 = z_o + 2*l*sin(a); 
	z_1 = z_o - 2*l*sin(a);
    /*----核心计算部分----*/	
	for(t=0; t<101; t++) {
		x_k[t] = t*20;	//Xk 0到2000 
		/*----参数部分------*/	
		r_1[t] = sqrt((double)(pow(x_k[t]-x_o+b+l*cos(a), 2.0) + pow(z_o-x_k[t]-l*sin(a), 2.0)));
		r_2[t] = sqrt((double)(pow(x_k[t]-x_o+b-l*cos(a), 2.0) + pow(z_o-x_k[t]+l*sin(a), 2.0)));
		r_3[t] = sqrt((double)(pow(x_k[t]-x_o-b+l*cos(a), 2.0) + pow(z_o-x_k[t]-l*sin(a), 2.0)));
		r_4[t] = sqrt((double)(pow(x_k[t]-x_o-b-l*cos(a), 2.0) + pow(z_o-x_k[t]+l*sin(a), 2.0)));
		
		v_1[t] = PI - Pan_d((z_o-z_k-l*sin(a)), (x_k[t]-x_o+b+l*cos(a)));
		v_2[t] = PI - Pan_d((z_o-z_k+l*sin(a)), (x_k[t]-x_o+b-l*cos(a)));
		v_3[t] = PI - Pan_d((z_o-z_k-l*sin(a)), (x_k[t]-x_o-b+l*cos(a)));
		v_4[t] = PI - Pan_d((z_o-z_k+l*sin(a)), (x_k[t]-x_o-b-l*cos(a)));
		/*------公式部分------*/
		san_g[t] = 2*G*Q*((double)((z_2*(v_2[t]-v_4[t])-z_1*(v_1[t]-v_3[t]))+x_k[t]*(sin(a)*sin(a)*log((r_2[t]*r_3[t])/(r_1[t]*r_4[t]))+cos(a)*sin(a)*(v_1[t]-v_2[t]-v_3[t]+v_4[t]))+2*b*(sin(a)*sin(a)*log(r_4[t]/r_3[t])+cos(a)*sin(a)*(v_3[t]-v_4[t]))));
		san_x[t] = (M/(2*PI))*sin(a)*((double)(log((r_2[t]*r_3[t])/(r_1[t]*r_4[t]))*cos(a-i) - sin(a-i)*(v_1[t]-v_2[t]-v_3[t]+v_4[t])));
		san_z[t] = (M/(2*PI))*sin(a)*((double)(sin(a-i)*log((r_2[t]*r_3[t])/(r_1[t]*r_4[t])) + cos(a-i)*(v_1[t]-v_2[t]-v_3[t]+v_4[t])));
		/*------写入文件------*/
		fprintf(fp, "%lf %lf %lf\n", san_g[t], san_x[t], san_z[t]);
	}
	/*------关闭文件------*/
	fclose(fp);
}

/*----关于atan(b/a)的判断,并返回一个角度------*/
double Pan_d(double b, double a)
{
	if(fabs(a)>1.0e-15) {
		if(b/a > 0.0) {
			if(a>0.0 && b>0.0) {
				return atan2(b,a);
			}
			else {
				return (-PI + atan2(b,a));
			}
		}
		else {
			return (PI + atan2(b,a));
		}
	}
	else {
		if(b<0.0) {
			return -PI/2;
		}
		else {
			return PI/2;
		}
	}
}